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Exploring Contractor Renormalization: Tests on the 2-D Heisenberg Antiferromagnet 

and Some New Perspectives * 

M. Stewart Siu^ and Marvin Weinstcin^ 
Stanford Linear Accelerator Center, Stanford University, Stanford, California 94309 

Contractor Renormalization (CORE) is a numerical renormalization method for Hamiltonian 
systems that has found applications in particle and condensed matter physics. There have been few 
studies, however, on further understanding of what exactly it does and its convergence properties. 
The current work has two main objectives. First, we wish to investigate the convergence of the cluster 
expansion for a two-dimensional Heisenberg Antiferromagnet (HAF). This is important because the 
linked cluster expansion used to evaluate this formula non-perturbatively is not controlled by a small 
parameter. Here we present a study of three different blocking schemes which reveals some surprises 
and in particular, leads us to suggest a scheme for defining successive terms in the cluster expansion. 
Our second goal is to present some new perspectives on CORE in light of recent developments to 
make it accessible to more researchers, including those in Quantum Information Science. We make 
some comparison to entanglement-based approaches and discuss how it may be possible to improve 
or generalize the method. 

PACS numbers: 75.40.Mg, 75.50.Ee, 02.70.-c, 03.67.Mn 
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1. INTRODUCTION AND OUTLINE 



Many problems in physics, in areas ranging from par- 
ticle and condensed matter physics to theoretical quan- 
tum computing, can only be treated by numerical meth- 
ods. Among them is the particularly interesting prob- 
lem of extracting the low energy behavior of a multi- 
dimensional system defined by a Hamiltonian with local 
interactions. While analytical methods can be applied to 
a few such Hamiltonians, existing methods generally re- 
quire enormous computational power to study systems of 
even modest size. For example, Quantum Monte Carlo 
can give highly accurate results for many systems, but 
its applicability can be limited by the fermion sign prob- 
lem, as well as the inaccuracies inherent in extrapolating 
finite size results to the limit of infinite volume. The 
most popular alternative, the Density Matrix Renormal- 
ization Group method (DMRG),lJ, is particularly suc- 
cessful in one dimension. With the advent of Quantum 
Information Science, it has been extended to simulate 
time evolution Q and multidimensional models (PEPS 
or tensor networks 0, Q). The idea underlying DMRG 
and its generalizations is a clever variational ansatz - the 
representation of states as contracted tensors. This ap- 
proach has many virtues. For example, it provides an 
upper bound on the ground state energy density. One 
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can also extract quantities such as total entropy which is 
approximately preserved in the subspace and use them 
to analyze phase transition 0, It has its limitations, 
however, as one has to deal with convergence issues and a 
very large number of states in more complicated systems. 
The closely-related method described in Q also has the 
requirement of a fixed finite lattice. 

In search for an alternative method, one might rea- 
sonably ask if there is a way that encodes the informa- 
tion about the ground state and low lying excited states 
not in states, but in operators . One method that does 
that is the Contractor Renormalization Group Method 
(CORE) [H, which, like DMRG, is an attempt to im- 
prove Wilson's real-space renormalizationfioj . Similar 
ideas have appeared in the past @, H| but we will follow 
the terminology of Q as it is the first general formulation 
of the method. (The same formulation was also indepen- 
dently proposed as Real-Space Renormalization Group 
with Effective Interaction (RSRG-EI) in [IH). While the 
method has intuitive appeal and has been used to study 
many collective phenomena (l2l. H3, H3| . unlike in the case 
of DMRG, we have relatively little grasp of exactly why 
and when the algorithm works. There seems to be a need 
to put aside applications for a moment and look at the 
method more closely. Short of very conclusive results, we 
present here some findings that may lead to a better un- 
derstanding of CORE. We approach the issue from two 
angles. First, we take the simplest 2-D model as a testing 
ground and carry out successive numerical renormaliza- 
tion using three different blocking schemes. Our choice, 
the Heisenberg Antiferromagnet (HAF), has been stud- 
ied with RSRG-EI in the past [111 , but our focus is on the 
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extraction of order parameter and the behavior of long- 
range operators in the 2-D cluster expansion, particularly 
because there is a freedom in defining successive terms in 
the expansion [l5j|. We also look at how well different 
blocking schemes agree with each other while capturing 
the physics in distinct ways. Secondly, we make some the- 
oretical and numerical comparison between CORE and 
DMRG as well as Entanglement Renormalization [lij ]. 
By presenting a slightly different perspective, we hope to 
provoke further investigation into the limitations of the 
method and the possibility for improvement. 

Section 2 of the paper reviews the basic formulas in 
CORE and sets the notation. The first subsection of Sec- 
tion 3 presents calculation of CORE on nine-site blocks, 
as considered in [ll|]. In addition to the energy density, 
which can be easily obtained to high accuracy, we calcu- 
late the staggered magnetization and find that, without 
longer-range operators, it is only accurate to one signifi- 
cant figure. To see how longer-range terms behave with 
limited computing resource, we use smaller blocks in sec- 
tion 3(b) (four-site) and 3(c) (five-site) and compute the 
energy density. We find that operators beyond nearest- 
neighbor can contribute quite significantly and requires 
very careful ordering. To this end we propose a ordering 
scheme based on diameters of the clusters. The effect of 
long-range operator, however, does not necessarily cor- 
respond to long-range correlation, as we see in 3(b) that 
even with only nearest neighbor operators a vanishing 
mass gap appears non-trivially. Although application is 
not the focus of the current work, we also show as a side 
note in Section 3(c) how, under the appropriate blocking, 
CORE might provide an interesting justification for the 
spin-wave approximation. 

In Section 4 we turn to the question of how CORE 
relates to entanglement-based approaches. In the first 
subsection we discuss truncation and blocking schemes 
of CORE and DMRG and compare their principles. We 
show that for a small, finite toy model CORE yields re- 
sults comparable to DMRG. Then in Section 4(b) we for- 
mulate CORE in a way that enables us to see its similar- 
ity to Entanglement Renormalization 16], which should 
be more familiar to readers in Quantum Information Sci- 
ence. We also discuss in Section 4(c) the role of block 
entropy in CORE and its possible use. Finally in Sec- 
tion 4(d) we discuss how the choice of retained states in 
CORE can affect its performance and what entanglement 
has to do with this choice. 

In Section 5 we reprise our results and discuss a number 
of future directions and possibilities. 

2. CORE: THE BASIC FORMULAS 

The original description of CORE can be found in 
Ref.@|. This section summarizes the basic formulas we 
will use in the sections to follow (particularly in Section 
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Given a Hamiltonian H on Hilbert space TL, the renor- 
malized Hamiltonian that CORE seeks to approximate 
takes the form: 

-^ren(0 — 

(p e -2Htpyi Pe -Ht He -m p ( Pe -2Ht p yi 

(1) 

where P is a projection operator from Ti onto W , a cho- 
sen subspace of retained states (in the language of RSRG- 
El[ll|, the model space), t is a variable parameter usu- 
ally taken to infinity. With the following lemma, we can 
show that H ren = lim^oo H ren (t) takes a particularly 
simple form. 

Lemma 1 Let |cij),i = l,...,M be an arbitrary or- 
thonormal basis for Ti.' and let \vi) ,1 = 1,...,N, be 
an orthonormal set of eigenvectors of the Hamiltonian, 
H , with eigenvalues {E\, E2, ■ ■ . , Ei . . .} arranged in as- 
cending order (generally, M < N). Then, there exists 
an M x M matrix, Rij, such that the states \wj) — 
^2iLi Rij \ a i) have the property that for each j , there is 
exactly one index 1 < f(j) < N such that (wj\vf(j^ ^ 
and (wk\vf(j-\) — OVfc > j . 

Corollary 1 Each state \wj) has the property that 

e~ m K) e"W (v m \ Wj ) \v f(j) ) + ..., (2) 

where . . . stands for terms that vanish more rapidly as 
t — ► 00. 

Then it can be proved that: 

= ^2\wi)E f(i) (Wi\. (3) 

i 

On a lattice we may expand this renormalized Hamil- 
tonian as linked clusters: 

H ren = Y, hConn ( s ^ 

sCL 

h conn (s) = H ren hC ° nn ( s ') (4) 

s'Cs 

where L is the entire lattice corresponding to Ti' (it can 
be infinite) and s is a connected sublattice. H ren \ s is 
H ren evaluated for the theory obtained by restricting the 
full Hamiltonian to the sublattice s. For notational con- 
venience, operators acting on s are implicitly assumed to 
extend to the full lattice by acting as the identity oper- 
ator on degrees of freedom lying outside of s. We will 
also refer to h conn (s) as a range-r connected operator if 
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s contains r blocks used to define the projection P, and 
we will adopt the notation h conn (s r ) to denote such an 
operator. 

In practice we can of course only evaluate range-r op- 
erators for small r, but since the longer-range terms are 
not controlled by a small parameter (this would be true 
if t is a small number) , we would have to be careful what 
to discard and we will discuss some examples in the next 
section. 



3. THE TWO-DIMENSIONAL HEISENBERG 
ANTIFERROMAGNET 

We have chosen the 2-D spin-1/2 HAF as testing 
ground for its simplicity and connection to the Hubbard 
model. The Hamiltonian is defined as 

H = ^2 Si - Sj 

<i,3> 

= X s * • s j ■ s t ■ s ': + s t ■ s h ( 5 ) 



<i,J> 



where the sum is over all neighboring pairs on a two- 
dimensional square lattice. 



a. Nine-site Square Blocking 

Although we are not able to compute (on a PC) long- 
range terms with nine-site blocks, their simplicity allows 
us to see some basic features of the computation and the 
HAF we are studying. The nine-site block is a natural 
extension of the three-site block in 1-D Q as it also has 
a spin-1/2 multiplet as ground states which we use as 
retained states. The range-1 and range-2 operators take 
the form: 



h conn (s l ) =a l 
l (s 2 ) = ail + 



(6) 



where we denote the one-block and two-adjacent-block 
configurations by s 1 and s 2 respectively. Note that no 
calculation of R in EqQ] is necessary to find h conn (s 2 ), 
because the retained states form exactly a spin-0 multi- 
plet and spin-1 multiplet. If we are to preserve the spin 
symmetry, there are no other inequivalent rotations. The 
fact that spin symmetry largely dictates possible eigen- 
vectors of the renormalized Hamiltonian greatly reduces 
computational effort, though arguably, we can only test 
the effect of CORE'S recipe for R in more complicated 
situations where there are many multiplets of the same 
spin. (From our experience it is also possible to break 
the spin symmetry and let CORE decide the best lin- 
ear combination. This often yields a good ground state 
energy but a poor spectral distribution.) 



With the range-1 and range-2 operators above, we 
can easily calculate the ground state energy density by 
summing a geometric series and reproduce the result of 
8, 11]. The energy per site obtained this way is —0.666 
and within 0.5% of the Monte Carlo result -0.669(l9j]. 
The situation becomes much more complicated, however, 
when we we proceed to calculate the staggered magne- 
tization, an order parameter of the HAF. There are two 
ways of calculating the expectation value of a renormal- 
ized operator. In [9] it is argued that other operators 
should take a form similar to Eq[3J 



O re 



R'OR, 



where now O is the matrix in the 



basis 



Oij = 



\o\ 



(7) 



(8) 



But since the staggered magnetization does not easily 
converge to a simple form as the Hamiltonian does, it is 
easier to calculate the expectation value by: 



(0) = 



8(H + JO) 
SJ 



1,7=0 



(9) 



. In other words, we add a multiple of the staggered 
magnetization operator, JM stagi to the Hamiltonian, use 
CORE to compute the ground state energy density, E(J), 
and extract the slope of this function at J = 0(see EqlH]). 
Once we have added JM stag to the starting Hamilto- 
nian the renormalized Hamiltonian is no longer a simple 
multiple of the original Hamiltonian and obtaining E(J) 
requires running the RG until it converges. In FigQ] we 
plot the ground state energy obtained in this way as a 
function of J. The staggered magnetization is the slope 
of the curve at J = and is obtained by fitting to a 
fourth-order polynomial in J. 

How small can we set the J-values to be? If we are 
allowed to use the knowledge that the energy obtained 
is accurate to two significant figures, we have to choose 
the J-values to be large enough so that the error in the 
slope will be at least one order of magnitude smaller than 
the slope itself. The expectation value thus obtained is 
M sta g = —0.29, whereas the exact Monte Carlo result 
is —0.30 [3]. We must caution though that if we use 
smaller J's or change the order of the polynomial fit, we 
can change the result by up to 20%. To calculate the 
order parameter more reliably, it seems that one has to 
find a way to capture more physical information in the 
renormalization. 



b. Four-site Square blocking 

In the case of the four-site block the lowest lying states 
form a spin-0 singlet and a spin-1 triplet. In this case, 
the renormalized Hamiltonian is no longer isomorphic, 
even at range-2, to the original Hamiltonian which had 
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FIG. 1: A plot of staggered magnetic field J vs energy density 
E. Staggered magnetization can be calculated from (M) = 
^^-jy—*- \j=o- The data points have to be close enough to 
see the curvature but sparse enough so that error in the energy 
will not affect the slope between the points too significantly. 

a single spin- 1/2 degree of freedom associated with each 
site. Furthermore, it appears to describe a theory with a 
non-vanishing mass gap, as subsequent RG-steps contin- 
ues to have a spin-0 and a spin-1 multiplet as its lowest 
energy eigenstates. This single-site gap, however, is just 
a reflection of the uncertainty principle cost one must 
pay for localizing the spin-1 excitation to a single block. 
By keeping these two multiplets at every step, we run 
the RG until the energy density converges and find that 
the gap between the two multiplets converges to zero. 
This reflects the fact that spin- 1/2 HAF is massless and 
agrees with the result of nine-site and five-site blockings, 
where the ground state is by contruction degenerate. We 
would like to note that this result is entirely non-trivial, 
since the original Wilsonian RSRG (the t = limit of 
the CORE formula for the renormalized Hamiltonian) 
predicts that the same theory dimerizes and has a non- 
vanishing mass gap. 

While this result gives us some confidence that even 
the lowest order cluster expansion is doing fairly well at 
extracting the correct physics, we obtain —0.710 for the 
ground state energy density - it is not nearly as accurate 
as the nine-site case. This is perhaps to be expected given 
that we have exactly diagonalized a mere 256 states here 
(whereas we have diagonalized 2 18 states in the nine-site 
calculation), yet in some sense it is still remarkable, since 
we have kept proportionally more states in every iteration 
(4-out-of-16 vs 2-out-of-512). To check the convergence 
of the cluster expansion, we add the operator correspond- 
ing to a sub-cluster consisting of four blocks arranged in a 
square. We do this only for the first RG step as the anal- 
ogous term in the next iteration requires diagonalization 
of 2 32 states. This already improves the energy density 
to —0.677 and allows us to obtain a staggered magneti- 




FIG. 2: Blocking of five-site stars. A total of five blocks are 
shown in the figure, range-2 terms are connected operators 
on two adajcent blocks, such as A and B. There are two types 
of range-3 terms: Examples of "L"s include ABC, BCD, etc., 
while there is a straight-line term on CDE. Finally, the range- 
4 plaquette acts on ABCD. 

zation that is close to the nine-site case. At this point we 
might be prompted to ask: To achieve a good accuracy, 
do we simply compute all the terms within our compu- 
tational power? Since the five-site blocking exhibits the 
most numerical sensitivity to the ordering in the clus- 
ter expansion, we defer a more detailed discussion to the 
next subsection. 

c. Five-site Star blocking 

The final CORE computation we wish to discuss uses 
the five-site blocks shown in Fig(2J What makes the RG 
transformation following from this blocking procedure in- 
teresting is that it behaves quite differently from the one 
obtained by restricting attention to square blocks. This 
is because the ground state of a star, in contrast to the 
nine-site block, is a spin-3/2 multiplet. The renormalized 
Hamiltonian at range two takes the general form 

h conn (s 2 ) = c l + Cl S ■ S + c 2 (S ■ S) 2 + c 3 (S ■ S) 3 . (10) 

If we evaluate the expectation value of EqlTO] in the clas- 
sical Neel state, we obtain an upper bound on the energy 
density of -0.712 per site. So despite the fact that the 
fundamental block has more sites, this accuracy is worse 
than the equivalent four-site blocking. Notice that the 
ratio of intra-block over inter-block links is smaller with 
five-site blocks. We suspect that the large perimeter of 
star blocks result in many under-constrained sites near 
the edge of the configuration. 
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Nevertheless, we can continue to run the renormaliza- 
tion group to see how the theory changes. If we use the 
Hamiltonian defined in Eq llOl and once again calculate 
the spectrum of a single star we find that now the ground 
state is a spin-9/2 multiplet. Keeping the spin-9/2 multi- 
plets as new degrees of freedom, we find that the range-2 
interaction again to be antiferromagnetic. We use these 
interactions to construct a five-site star, and the ground 
state now becomes a spin-27/2 multiplet. From this ver- 
sion of CORE we see that the spin-1/2 theory is equiva- 
lent to a theory with a larger spin at each site. 

Clearly the stability of a picture that says that spins 
keep growing with each renormalization group step must 
be checked by computing the contributions to the renor- 
malized Hamiltonian coming from larger clusters. The 
reason this is interesting is that once one goes to larger 
sub-clusters, one introduces new diagonal couplings, cou- 
plings involving three sites at a time, etc. Since the diag- 
onal couplings are also mainly antiferromagnetic in char- 
acter, they tend to introduce frustration into the system. 
It is entirely possible that these are relevant operators 
and significantly modify the nature of the RG flow, so 
that at the next step the spins fail to grow. 

To obtain a qualitative understanding of how each con- 
nected operator modifies the flow of the renormalized 
Hamiltonian we first observe that each three-site con- 
nected operator can be written as a sum of terms which 
act non-trivially on one, two or three-sites at a time. To 
see why this is true let us assume that we have spin-n/2 
on each site of the new lattice and let Xi,i = 2..n + 1 
be a set of n traceless matrices which, together with the 
matrix X\ — l/(n + 1) 1 (where 1 is the unit matrix), 
form a basis for the space of of (n + 1) x (n + 1) matri- 
ces. Furthermore, assume that these matrices satisfy the 
normalization conditions 

TH.V,.\j d,,. (11) 

It then follows that the tensor products 

M ijk =Xi® Xj (g) Xj (12) 

form a normalized basis of matrices which operate on 
three lattice sites at a time. Thus, the general three-site 
connected operator can be written as 

O = Ys a vkMi jk . (13) 

Furthermore, the coefficents aijk can be extracted by tak- 
ing the trace with any M^k] i.e., 

Oi jk = \ v:.\l, ,,(),. (14) 

Given this definition, it is easy to define what we mean 
by the parts of O which act non-trivially on zero, one, 
two or all three sites. The part of O that is proportional 
to a multiple of the unit matrix is, of course, a term that 



acts trivially on all three sites and won't contribute to 
dynamics. Of what remains, it typically turns out that 
the operators which act non-trivially on only one or two 
sites are the most important operators. Thus, for the 
purpose of getting a simple qualitative understanding of 
the stability of our problem we restrict attention to these 
operators in what follows. Since the renormalized Hamil- 
tonian must commute with the total spin operator, it 
follows from direct computation, or simple group theory, 
that there are no terms in the renormalized Hamiltonian 
which act non-trivially on a single block. (This is because 
for spin-n/2, the space of (2n + 1) x (2n + 1) matrices 
decomposes, under total spin, into matrices which trans- 
form as spin-0, the unit matrix, and, for each j, a set of 
matrices which transform as spin-j, where j runs from 
1...2n+ 1.) Similar arguments tell us that operators 
that act non-trivially on only two-sites at a time can be 
written as polynomials in S ■ S acting on the two-sites 
in question. By S we mean the matrices which represent 
the generators of spin rotations for spin-n/2. Note, the 
order of the polynomial for the case of spin-n/2 is at most 
2n. 



Having said this we can ask what the effects of the 
terms proportional S ■ S (which act on only two-sites at 
a time) coming from the straight line, the L and the pla- 
quette are (this terminology is best explained by looking 
at Fig[2]), since these two-site operators turn out to be 
a significant part of the contribution to the renormal- 
ized Hamiltonian. There are three important observa- 
tions that we must make about these terms. 



First, the antiferromagnetic interaction between di- 
agonal sites (e.g., BD in Figj2|) are of the same or- 
der of magnitude as the interaction between adjacent 
sites. This means that the contribution of larger clusters 
to the renormalization group transformation gcncrically 
produce a frustrated system. Second, unlike the case of 
one dimension, where the importance of operators com- 
ing from larger clusters falls off with the number of blocks 
in a cluster, terms coming from the four-block plaquette 
appear to be as important as those coming from clusters 
consisting of only three-blocks. (Among the three types 
of terms we evaluated, the three-block straight-line has 
the smallest two-site contributions.) Third, the two-site 
terms coming from different clusters often cancel each 
other, indicating that our result is sensitive to how we 
define the cluster expansion. To give the reader a feeling 
for these three points we include the following table which 
shows how the energy density changes with the inclusion 
of different terms when we approximate the ground state 
with the Neel state at the spin-3/2 level. 
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Terms included 


Approx. energy density 


Only range-2 


-0.712 


+"L" 


-0.623 


+"L"+"line" 


-0.586 


+" L" +plaquette 


-0.637 


+" L" +" line" +plaqucttc 


-0.600 



At this point it is natural to ask, "What is the best 
way to arrange the terms in the cluster expansion?" . In- 
tuitively, we expect the contribution of a long chain in a 
straight line to be less than the contribution of a square 
with the same number of blocks, because the blocks in the 
latter case are closer to each other. Taking this with the 
observations above, we propose a cluster re-summation 
that treats all clusters having the same diameter as hav- 
ing equal weight. By the diameter of a cluster we mean 
the longest distance between two blocks in the configu- 
ration. For example, the diameter for the range-2 term 
would be one, for the three-block "L" and four-block pla- 
quette term y/2, etc. This scheme has the additional 
advantage that now each term in the cluster expansion 
preserves at least some of the rotational symmetry of the 
lattice since, at every order, we include all the terms with 
the same diameter. 

Note that in a diameter expansion the four-block pla- 
quette comes before the three-block straight-line, which 
has diameter two. Of course, to be sure that all terms 
coming from clusters of diameter 2 are smaller than those 
of diameter 1 and diameter y/2, we should have also eval- 
uated the cluster where four blocks are arranged in a 
"T". We left it out because, for symmetry reasons, this 
configuration is considerably more difficult to evaluate 
than the plaquette and we wanted to limit this study to 
computations easily carried out on an average personal 
computer. Using only the plaquette and "L" terms, the 
ground state of the five-site spin-3/2 star again turns out 
to be a spin-9/2 multiplet. This is somewhat remarkable. 
Even with pure S ■ S type interactions between adjacent 
and diagonal sites, the ground state of the generically 
frustrated five-site star could have been spin-1/2, spin- 
3/2 or spin-9/2 depending on relative strengths of the 
interactions. 

In principle we should now proceed to calculate the 
interactions for the spin-9/2 theory up to the same di- 
ameter. Unfortunately, calculations of the plaquette in 
the spin-9/2 theory requires the diagonalization of a 
sparse 4 20 x 4 20 matrix, so we are restricted to range-2 
(diameter-1) as in the four-site case. The range-2 inter- 
action remains antiferromagnetic but by itself gives an 
energy density that is less than —0.78. 

Although, due to limited resources, we do not have 
the longer range terms to show a converging result, we 
would like to close this discussion by noting that this 
is very interesting in the context of the spin-wave ap- 
proximation to the HAF. Since turning the spins into 
Holstein-Primakoff bosons rely on an expansion in 1/2S, 
a validated picture of growing spin essentially explains 



the well-known accuracy of the spin-wave method on 
the spin-1/2 lattice. We can estimate the error in the 
1/2S expansion for a finite spin-9/2 lattice. We take 
our range-2-only calculation and rewrite the Hamilto- 
nian on one five-site block of spin-9/2 using Holstein- 
Primakoff bosons, keeping only oscillator terms that are 
of no higher power than the number operator. We then 
find the ground state energy by minimizing over Bogoli- 
ubov transformation parameters. The value calculated 
this way differs only by 0.02% from the exact energy of 
the spin- 2 7/ 2 multiplet. 

In conclusion, the three-blocking schemes give us three 
different physical pictures for the same system. The sen- 
sitivity of CORE to link structures suggests that, when 
working with exotic blocking schemes, one has to check 
very carefully the convergence of the cluster expansion. 
Moreover, our experimentation with various ways of re- 
summing range-3 and range-4 terms indicates that the 
RG-flow can be changed dramatically depending on how 
we group and order the operators in the cluster expan- 
sion. The diameter expansion we proposed appears to 
be the most plausible solution, but in absence of rigorous 
theorems bounding the long-range terms, this remains an 
open problem. 



4. COMPARISON TO ENTANGLEMENT-BASED 
APPROACHES 

a. Density Matrix Renormalization 

The remainder of this paper is devoted to comparing 
CORE to other methods in a way we hope would shed 
light on the features of CORE. Given the popularity of 
DMRG, we take it as an instructive benchmark for study- 
ing CORE's performance. This sort of comparison has 
not been made in the past because the two methods in 
their respective original formulations have very different 
blocking and truncation schemes and it seemed difficult 
to compare them in a meaningful way. There are two 
essential features in the original formulation of DMRG 
(we refer the readers to Refs.[l[ for details): 

I. Reduced Density Matrix Truncation - In each 
block we select what to keep according to the 
reduced density matrix of a target state, usually 
the ground state of a larger system. The larger 
system is the block whose truncation we are 
interested in plus an " environment" , which can 
be for example a copy of the block. Truncation 
consists of eliminating those vectors with the 
smallest eigenvalues of the reduced density matrix 
obtained by tracing out the environment. The 
error in the truncation depends on the distribution 
of eigenvalues of the reduced density matrix and 
can be analyzed in the language of quantum data 
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FIG. 3: A) Hierarchical Blocking allows renormalization in 
the Wilsonian sense. B) Linear Blocking instead grows the 
block by one site at a time and works better for the original 
DMRG. (Details about the environment omitted above.) 



compression [17j • 



II. Linear Blocking - The block size increases linearly 
one site at a time. This naturally allows iterative 
improvement by back and forth sweeps and gives 
rise to the underlying matrix product state struc- 
ture. Yet this is not renormalization in the sense 
of coarse-graining one theory to another of a differ- 
ent scale. To do so one could, of course, use hier- 
archical blocking (Fig[3J\.) along with the reduced 
density matrix truncation described in I (We will 
refer to this as "Hierarchical DMRG" or HDMRG), 
but it is usually not preferred. Roughly speaking, 
the reason is that entanglement often scales with 
the boundary of the block, and because hierarchical 
blocking gives rise to more boundaries, truncation 
in HDMRG results in more loss of information. 

CORE as formulated in Section 2 does not rely on spe- 
cific truncation and blocking schemes - they are details of 
the projection operator P. The scheme we use in Section 
3 can be classified as "hierarchical blocking" (FigJHJ'V) and 
we will refer to it as such. (Though unlike HDMRG or 
naive real-space renormalization, the state cannot be re- 
constructed in the original space in a simple hierarchical 
manner.) 

Since we are free are to choose the form of the projec- 
tion operator, this difference between CORE and DMRG 
can be considered a superficial one. As we mentioned in 
the Introduction, the essential feature of CORE is that it 
encodes some information in operators instead of states. 
To see this in a formulation that allows direct compari- 
son, let us consider the alternative form of DMRG - the 
matrix product states (MPS) formalism. One common 
way to obtain the ground state energy using MPS is by 
applying an imaginary time evolution to a simple starting 
state, which we will call |^). Suzuki- Trotter decomposi- 



tion can be used to decompose the evolution into small 
steps and after every step the state will be truncated to 
make sure that it lies within the subspace spanned by 
MPS of a fixed dimension. If the procedure converges 
successfully to the desired attractor and P is a projec- 
tion onto the MPS subspace, the final state should be 
the same as Pe~ Ht \tp) = Pe~ Ht P \ip) where t is taken 
to infinity and the equality follows from the fact that the 
starting state lies within the projected subspace. The 
ground state energy is therefore: 



En 



Jj\Pe- m PHPe- Ht P\ij) 
(iP\Pe- Ht Pe- Ht P\ip) 



(15) 



This means Eq is the smallest eigenvalue of the Hamilto- 
nian: 



(Pe- m Pe- nt P) 



-Ht - 



-?Pe- Ht PHPe- Ht P(Pe- Ht Pe- Ht Py 



(16) 



We can now contrast this directly with Eq[TJ H ren (t) 
in Eq[T]contains the exact ground state energy but cannot 
be evaluated exactly, therefore we have to throw away 
the long-range terms caused by the diffusion of e~ Ht . 
When the additional projections are inserted in H' ren (t) 
of Eq[l6l it no longer contains the exact ground state 
energy, but its ground state energy can be found exactly. 
In this case the burden of a good approximation is shifted 
entirely to a clever choice of P. 

Having seen an abstract comparison, let us now re- 
turn to the original form of DMRG and compare some 
numbers. Because hierarchical blocking allows CORE 
to handle very large lattices and DMRG requires linear 
blocking, one might expect CORE's strength to be with 
large systems. We were surprised to find, however, that 
even on small finite lattices, CORE achieves accuracies 
that are comparable to DMRG. We demonstrate this by 
running HDMRG, DMRG and CORE on a short periodic 
Ising chain. 

The first set of plots, FigfJ] is for the nine-site peri- 
odic Ising chain and compares CORE to HDMRG. This 
is an interesting comparison since we can naturally use 
hierarchical blocking for both. The first plot is for the 
ground state energy and the second one is for the gap be- 
tween the ground state and the first excited state. The 
periodic chain is divided into three blocks and we retain 
two states per block. To produce the HDMRG plot, we 
exactly diagonalize all nine-sites to generate the target 
state. In other words, we take the exact ground state, 
trace out six sites and use the reduced density matrix on 
three sites to determine what to keep in each block. At 
the end we diagonalize the Hamiltonian truncated to the 
retained states to produce the numbers for the plot. Of 
course, in a realistic calculation, we do not have the exact 
ground state to use as target state, but it can be approx- 
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imated by iteration, so we are essentially simulating the 
ideal converged limit. 



CORE calculations are carried out up to range-2 in 
the same manner as in Section 3. (Since we are deal- 
ing with a finite lattice, range-3 would be exact.) In 
the renormalized space we'd have three sites and a new 
nearest-neighbor Hamiltonian, which we can use to find 
the energy and the gap. Actually, there are two versions 
of CORE in our plots with different P operators. We 
will elaborate on this further in Section 4(d). The first 
version (represented by triangles) retains the two states 
with the lowest energies in each block just as in Section 
3. The second version (represented by circles) uses a 
reduced density matrix truncation scheme: it takes the 
ground state of six sites, trace out three sites on one side, 
and retain the two states with the highest weight in the 
reduced density matrix. The plot indicates that the sec- 
ond version performs a little better in energy but at the 
expense of the gap. Note that the six-site calculation is 
only one out of many ways to obtain a reduced density 
matrix. In fact if we use the ground state of all nine sites 
as in our HDMRG calculation instead, the states with 
the highest density matrix weight turns out to be the 
same as those of the first version. We have chosen the 
six-site case because they are what we use for the range-2 
operator calculation. 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Explanation of Symbols/Corresponding Blocking 



Triangles: Range-2 CORE 
with lowest energy truncation, 
2 states/block 

Circles: Range-2 CORE with 
reduced density matrix 
truncation, 2 states/block 
Crosses: "Hierarchical 
DMRG" 2 states/block 




FIG. 4: Performance of various methods on a nine-site peri- 
odic Ising system. The two CORE cases (triangles and circles) 
use different truncation schemes. One keeps the lowest energy 
states within the block; another finds the ground state of two 
blocks and considers its reduced density matrix. They are 
compared to HDMRG with the same blocking scheme. 



FigfJ] shows that CORE compares favorably against 
HDMRG, so let us next compare CORE to DMRG with 
a blocking scheme that is natural to the latter. The cho- 
sen system is an eight-site periodic Ising chain with two 
blocks and two free sites, shown in Figj5] The two three- 
site blocks mimics the system/ environment in the middle 
of a DMRG sweep, and the two free sites are positioned 
according to the prescription in Ref . [23| . Two states are 
kept in each block. As in the HDMRG case, we simulate 
the iteration (again we refer the readers to Refs.[l|, [23| 
for details) by using the exact ground state as our target 
state, so our result represents the convergence limit. If we 
had actually carried out the iterative sweeps, we would 
need to diagonalize at most 2 4 states at a time (two states 
in each free sites, two states in the system block and two 
states in the environment block), so for fairness, we cal- 
culate CORE only up to range-2 because it also requires 
diagonalization of 2 4 states (one block plus one free site 
exactly). Here we simply use the first version of CORE 
which truncates to the lowest energy in each block. Plots 
in Fig[5]show that both DMRG and CORE achieve very 
high accuracies, but the two methods excel at different 
regimes of the Ising system. This indicates that CORE's 
performance, at least in simple situations, is comparable 
to DMRG without any need for iteration. 
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Explanation of Symbols/Corresponding Blocking 



Pluses: Range-2 CORE, 
2 states/block 
Diamonds: DMRG, 
2 states/block 




FIG. 5: Comparison of CORE to DMRG on an eight-site 
periodic Ising system. The blocking mimics the intermedi- 
ate configuration of a DMRG sweep with periodic boundary 
condition. 



The presence of entanglement is often believed to be 
what makes quantum systems difficult to simulate. By 
"entanglement" we mean a measure that quantifies how 
much a quantum state deviates from a tensor product 
of states on subsystems. In the case of a pure state, 
for example, we can measure how much a finite block is 
entangled to the rest of the system using its von Neumann 
entropy, 



S(p) = -Tr(plnp) 



(17) 



where p is the reduced density matrix of the block. Effi- 
cient preservation of entanglement accounts for much of 
DMRG's success [ij. 

Without explicit consideration of entanglement, how 
does CORE achieve a performance comparable to 
DMRG? In this subsection we will try to understand this 
through theoretical comparison with an entanglement- 
based method. 

Vidal recently proposed a method called Entanglement 
Renormalization (ER) [l6| , which we can think of as an 
improvement on HDMRG. Recall that, by HDMRG, we 
refer to the method that uses reduced density matrix 
truncation as DMRG does but blocks hierarchically as 
CORE and naive renormalization do. Vidal observed 
that before we truncate states in a block according to 
the reduced density matrix of a target state, it is possi- 
ble to apply a unitary transformation on the boundary 
sites of two adjacent blocks and reduce the entanglement 
of the target state. In other words, the transformation 
makes states which are truncated away have less weight in 
the reduced density matrix of the block. Thus the same 
number of retained states can preserve more of the infor- 
mation of the original system. In exchange, we pay the 
price that operators acting on one block act on neighbor- 
ing blocks after the disentangling unitary transformation. 
The renormalized Hamiltonian can be written as: 



H?£ = W^UHU^W, 
U = U X .2 ® U 2 , 3 <8> ... 



(18) 



where the {t/^ i+1 }'s are disentangling unitary transfor- 
mations acting on the edges of neighboring blocks + l 
and W is an isometric transformation (W-^VF = /) lift- 
ing from the coarse-grained subspace to the full Hilbert 
space. For example, if {|cii)} form a basis of the sub- 
space and {\bi}} form a basis of the full space, we can 
write Wij = (bi\a,j), i.e. a matrix whose columns are the 
basis vectors of the subspace. The orthogonal projection 
P is related to W by P = J2 k \ a k) (a k \ = WW^ . 

Like CORE, Entanglement Renormalization stores 
some information in operators in addition to the retained 
states. The two methods turn out to be have a similar 
starting point. In Section 2 we have shown the form of the 
renormalized Hamiltonian in the limit of t — > oo. Now, 
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since both and (denned in LemmaUJ are 

orthonormal, we can think of the mapping which identi- 
fies each {|wj)} with its corresponding {^/(j))} as the 
restriction of a unitary transformation which acts on the 
full Hilbert space. Of course, given our construction we 
only have the restriction of this transformation to the 
subspace spanned by the retained states, the extension 
of the transformation to the full Hilbert space remains 
undefined and is presumably not unique. The fact that 
this unitary transformation is not uniquely specified is 
equivalent to the fact that W and W\ which appear in 
the ER formulas, are isometries and not unitary trans- 
formations. 

To put the correspondence between CORE and ER in 
a formal setting we make the following claim: 

Claim 2 For all unitary operators U such that 

u \ v f(j)) = K)> 

H ren = PUHU^P. (19) 

Proof: 

pu = 53 K) (vh\u 

i 

i,l 

= ^2\ w i) s f(i)i( v i\ 

i,l 
i 

(20) 

where the third line follows from the observation that 
U \ vi) must be orthogonal to Ti! if Vi, I ^ /(*)■ Hence 
we can write 

PUHtfP = J2\ w *H v f(i) \ H \ v fU)) ( w j\ 

= Y^\wi)E m (wi\ (21) 

i 

which is exactly Eq[3j ■ 

As we noted earlier, we can write P = WW' . In prac- 
tice, we want to write H ren in a basis of the subspace 7i', 
so what we really calculate is H ren = W'UHWW , just 
as in EqfTBJ Thus, we see that both CORE and Entangle- 
ment Renormalization use a renormalized Hamiltonian of 
the form PUHU^P. The distinction between CORE and 
ER is that Entanglement Renormalization approximately 
disentangles a system, while CORE approximates a dis- 
entangled system. We say that Entanglement Renormal- 
ization approximately disentangles the system because 
there is no guarantee that we can find disentangling uni- 
taries that reduce the rank of the reduced density matrix 



to less than or equal to the dimension of retained states. 
It is usually necessary to truncate some states to keep 
the degrees of freedom manageable, and information is 
lost when we truncate the states. CORE approximates a 
disentangled system in the sense that we first write down 
the form of a completely disentangled system fEq [T9")) . Its 
cluster expansion truncated to diameter-fc then approxi- 
mates this system by ensuring that the new Hamiltonian 
restricted to any sublattice with diameter less than k is 
completely disentangled. Here information is lost when 
we truncate the operators with diameter larger than k. 

In this picture, the CORE prescription for choosing 
R and Ef(j\ is a particular choice of disentangler. One 
might conceive of a disentangler PU that does not always 
require overlaps between \wi) and - a trivial exam- 

ple would be to set H ren = £\ \ a i) E ( a i\ where { |a^) } 
is the basis of H' that we start with. The problem with 
such an arbitrary choice is that the exact renormalized 
Hamiltonian could be so non-local that it is difficult to 
approximate it by a cluster expansion. Our intuition is 
that by keeping \wi) somewhat "close" to cluster 
expansion can do well. There may be, however, room for 
improvement once we have a better understanding of the 
cluster expansion. 

c. The use of entropy in CORE 

The von Neumann entropy (EqfTT]) of finite blocks in 
a lattice is known to exhibit scaling behaviors with block 
size that depend on the dimension and phase of the sys- 
tem Q. Because DMRG truncates according to the re- 
duced density matrix, entropy can be approximately pre- 
served, so apart from using entropy as a measure of how 
much information is lost, it is also possible to use the ap- 
proximate entropy measures from DMRG to detect phase 
transitions [||. 

Does entropy have a similar use in CORE? Note that 
CORE is not a variational method that works with a 
subspace within the original space. It tries to preserve 
eigenvalues but not the eigenvectors and this allows a 
disentangling effect. Therefore there is no a priori reason 
to expect CORE to preserve any particular amount of 
entropy 

Nonetheless, when we evaluate the three-site block en- 
tropy for the two toy models in Section 4(a), we were im- 
pressed by how much information CORE captures. Fig[6] 
shows the entropy of some blocks within the eight-site 
and nine-site models (the symbols are explained Fig[4] 
and Figf5]). The first two plots are for a three-site block, 
in which the eight states are truncated to two after the 
renormalization. This means that the upper bound of the 
block entropy changes from In 8 to In 2. As it turns out 
the exact block entropy in the eight-site case is far from 
saturating the bound, allowing DMRG and CORE to ap- 
proximate it closely. The third plot in FigJS] is for the 
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Von Neumann Entropy of Three-site Block for the Nine-site Ising Model 
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FIG. 6: Here we show the block entropy as a function of cou- 
pling for the nine-site and eight-site Ising models. The solid 
curve indicates the exact entropy and other symbols represent 
methods explained in previous plots. The first two plots are 
for three-site blocks and the last for a single-site. We see that 
CORE reproduces the shape in a manner similar to DMRG. 



single-site entropy which is available only for the eight- 
site configuration (The sum of single-site entropies can 
be used to define the total information encoded in the 
wavefunction [ItI] ) . Here there is no truncation in the site 
itself, so the entropy upper bound is the same before and 
after the renormalization. 

It may seem confusing how CORE can have a dis- 
entangling effect and yet keeps a block entropy that's 
comparable to DMRG. By "disentanglement" we mean 
a reduction of weight on the eigenvalues of the density 
matrix outside of the retained space, which does not say 
anything about the distribution of eigenvalues inside the 
retained space. It is the latter that determines the block 
entropy after truncation. 

What is remarkable about FigJSl however, is not the 
amount of entropy CORE keeps, but the fact that it 
varies smoothly with the coupling A and has a shape sim- 
ilar to the exact entropy. This raises the possibility that 
entropy in CORE can also detect quantum phase tran- 
sitions. Regrettably, unlike Entanglement Renormaliza- 
tion we cannot perform disentangling unitaries without 
truncation to check entropy scaling in critical and non- 
critical systems 16], which means our result can be highly 
dependent on our choice of retained states and cluster ex- 
pansion. In absence of more data, we are hesitant to draw 
conclusions at this point, but this can be an interesting 
area to explore. 



d. The choice of retained states 

In Fig|4l we have shown two versions of CORE with 
different choices of retained states, one of which is in- 
spired by DMRG. Although we have not done so, for 
the study of entropy we just discussed one can even try 
other variants of DMRG that choose retained state s by 
maximizing entropy [20] or minimizing the Holevo x [l7j |. 
How choices of retained states affects the performance of 
CORE have never been investigated in the past so we 
would like to briefly discuss it here. 

In the single-block calculation it is clear that whatever 
states we retain, as long as they have an overlap with the 
lowest lying states of the block, the renormalized range-1 
(i.e. single site) Hamiltonian will be unchanged. It is only 
when we compute the connected range-2,range-3, etc., 
operators that the difference in choice of retained states 
show up. In general, the effect of changing the choice of 
single block retained states can either increase or decrease 
the size of the longer range connected interactions. 

Past works on CORE have exclusively retained states 
with the smallest energy in the single block Hamiltonian. 
Let us call this the first version of CORE. Inspired by 
DMRG, we have tried a second version that uses states 
with the largest eigenvalues in the reduced density matrix 
of some target state. When we choose the target state 
to be the exact ground state of two adjacent blocks, we 
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often obtain a slightly better ground state energy at the 
expense of a less accurate gap between the ground state 
and the first excited state, as can be seen in FigJU When 
we choose the target state to be the exact ground state 
of a larger block that contains our block in the center, 
however, the retained states often turn out to be the 
same as those in the first version, i.e., the states with 
the smallest single block energy. This latter observation 
is in agreement with the results of Capponi et al [3], 
who found that the retained states in the first version 
have very high weights in the reduced density matrix of 
a target state of a larger block. 



Capponi et al considered this question in a different 
context and proposed that one should use the reduced 
density matrix to check if the truncation in the first ver- 
sion of CORE is justified. We think that such a diag- 
nostic tool has to be used very cautiously. To see why 
we say this, consider first the case when the target state 
is calculated on two adjacent blocks. Here the diagnos- 
tic tool essentially measures the overlap between 
and \wi) (on the Hilbert space of two blocks) in EqfT9l 
A large overlap, however, does not guarantee that U is 
close to the identity, let alone a better convergence of 
the cluster expansion. An alternative is to use the mixed 
state corresponding to combination of all as the 

target state, which should more meaningfully measures 
how much U differs from the identity. Yet it still does 
not in general reliably indicate how well the cluster ex- 
pansion converges. All we can be certain of is that in the 
limit U —> I, the cluster expansion tends to be exact at 
nearest neighbor range; this does not tell us much when 
U is significantly different from /. In particular, if we 
consider the case when the target state is not the ground 
state on two blocks, but on a larger block containing the 
block we want to truncate, the situation is even less clear 
because in that case CORE does not explicitly rotate this 
ground state to the retained states. 



The reduced density matrix tells us a lot in DMRG, but 
as we have mentioned in the previous subsections, CORE 
is performing some disentangling action and thus does 
not preserve the entanglement structure of the states as 
DMRG does. Simply preserving the most entanglement 
in the original Hilbert space does not guarantee the most 
accuracy in CORE. Fortunately, it turns out that in our 
Ising example, when we choose our retained states to 
maximize (^/(i) l^i) on two blocks, we do get a better 
overall ground state (and sacrifice the accuracies of the 
excited states). This seems to indicate that the reduced 
density matrix can be helpful - as long as we remember 
that even as (vf^wi) goes to 1, there is no substitute 
for the longer-range terms in the cluster expansion. 



6. DISCUSSION AND FURTHER DIRECTIONS 

This paper has presented a number of numerical results 
on various models. While they have shown an usefulness 
consistent with the considerable body of literature on 
CORE/RSRG-EI, the picture is still far from clear. The 
fact that we get an improvement in energy from longc- 
range operators in 2-D shows non-trivially that there is 
a genuine small parameter hidden in the cluster expan- 
sion, which we believe is associated with the diameter 
and not the number of blocks in the configuration. Yet 
we still have little knowledge how this convergence man- 
ifests itself under different blocking schemes. When we 
tested the five-site blocking shown in Section 3(c), we 
expected the growing spin picture to be valid on prior 
physical grounds, so we were surprised by how much the 
long-range terms could affect the result. There seems 
to be a need for a way to at least estimate the effect of 
long-range operators. As far as we know, only a very few 
applications of CORE/RSRG-EI have considered long- 
range operators up to " diameter-\/2" [13j. This is pre- 
sumably due to limited computational power. Is it possi- 
ble to use other approximate methods of diagonalization, 
such as perturbation theory or DMRG, to estimate these 
long-range operators? This hybrid approach is certainly 
a direction we would like to pursue in the future. 

The second part of the paper, where we compare 
CORE to entanglement-based approaches, also raises 
a number of questions. In Section 4(b), we showed 
that CORE is in fact theoretically similar to Entangle- 
ment Renormalization. Naturally it would be interest- 
ing to see how ER performs numerically in comparison 
to CORE, and in particular, whether one can reproduce 
the three pictures of the antiferromagnet with ER. Sec- 
tion 4(c) raises the possibility that apart from the energy 
spectrum, we may be able to use the entropy to study 
phase transitions with CORE. This, if true, would be 
remarkable given that CORE was not designed to pre- 
serve entropy. In Section 4(d) we mentioned the use of 
reduced density matrix as an alternative method of decid- 
ing what states to keep in each block. There is a special 
circumstance when it can be important. This has to do 
with a situation which comes up as soon as one stud- 
ies frustrated HAF's where, with increasing frustration, 
single block levels tend to cross. Since one expects that 
the most interesting physics of these models is associated 
with these regions one has to know what states to choose 
when degeneracies occur. The most obvious solution is 
to keep all states which can cross as a function of the 
truncation, however this will increase the complexity of 
the numerical computations. In that case it may be ad- 
vantageous to retain the states determined by a DMRG 
calculation which have the highest weight and correct 
spin. 

Finally, it would be remiss of us not to point out, with- 
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out giving details, the possibility of adapting the princi- 
ples of CORE to other applications. We noted in Section 
4 that specific truncation and blocking methods are de- 
tails of the projection operator; the underlying principle 
is a general one about simplifying states at the expense of 
generating non-local operators. The principles of DMRG 
have been generalized to simulate real time evolution [2| - 
Can the principles of CORE can be applied to time evo- 
lution as well? There are good motivations for thinking 
about this. One of the most efficient simulation method 
for a specific class of unitary evolutions is the stabilizer 
formalism 25j , where we do not keep explicit information 
about the states and instead keep track of them using a 
set of operators. Since CORE allows one to calculate ex- 
pectation values at the expense of state information, we 
could ask if cluster expansion can be similarly efficient 
for certain types of unitary evolution. Apart from a di- 
rect simulation of unitary evolution, we may also consider 
turning a unitary evolution problem into a ground state 
problem using results in quantum complexity theory [24j| 
(Linear, instead of hierarchical, blocking would have to 
used in that case). All such possibilities would be inter- 
esting to explore. 
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Appendix: QR-Decomposition and Proof of 
Lemma [TJ 

While constructive proof of Lemma [T] has been given 
in the past 9], for practical implementation, we have 
found recursive QR-Decomposition to be a fast and con- 
venient way of calculating the rotation R (particularly 
when packaged libraries such as NAG or LAPACK are 
available). Given any M xN matrix A, M < N, the QR- 
Decomposition is defined by QTZ = A, Q T Q =1 where Q 
is an M x M matrix and TZ is an upper-triangular M x N 
matrix. For our application, A is simply the overlap ma- 
trix On , and the matrix Q plays the role of W (Thus the 
TZ of QTZ is not the R of Lemma [TJ. Note that the upper 
triangular matrix does not necessarily satisfy the contrac- 
tion condition (Eq[2j; it merely guarantees zeros entries 
below the diagonal {TZ(i,i) | i = 1..M}, a particularly 
weak condition if M <C N. Our solution is to start from 
the upper-left corner, move down the row and column 
after every QR-Decomposition until we find another non- 
zero entry that has non-zero entries below it (this means 
more than one retained states contract to the eigenstate), 
at which we perform another QR-Decomposition on the 
submatrix with that entry as the upper-left corner. We 
repeat this recursively until the submatrix size is zero. 
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